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ABSTRACT 


Optimal age replacement policies are designed to cut down system failures and 
minimize maintenance cost. By scheduling planned replacements, a system is replaced 
at age ¢ or at failure, whichever comes first, and the cost of replacement before failure 
(planned) is less than the cost after failure (unplanned). In this thesis, the distribution of 
lifetimes is a known, increasing failure rate phase type distribution. To find the optimal 
age of replacement, the parameters of the underlying phase type distribution must be 
estimated. 

An optimal age sequential estimation procedure is developed. In particular, the 
phase type distributions parameters are estimated using a matching moments nonlinear 
programming approach. Since there are many parameters associated with phase type 
distributions and the distributions include matrix exponential terms, the parameters are 
in general difficult to estimate. A specific case where the phase type distribution has 
initial probability vector ~a=(1,0,0) is studied for different sample sizes and compared 


with a similar nonparametric procedure. 
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I, INTRODUCTION 

Normally a system or component is replaced when it fails and C,, the cost of 
replacement before failure (planned replacement), is often less than C,, the cost after 
failure (unplanned replacement). The problem of determining when to repair and when 
to replace failing systems is the concern of management resources. Inefficient 
management due to the use of nonoptimal maintenance policies can lead to significant 
system maintenance cost. In general, optimal maintenance policies are designed to cut 
down the number of system failures and minimize maintenance costs. In this thesis, we 
study the problem of estimating a type of optimal maintenance policy (the optimal age 
replacement policy) when the underlying system life distribution is phase-type. In 
particular, we consider an adaptive estimation scheme where the estimated optimal 
replacement age is updated each time the system is replaced. 

Maintenance is defined to be all activities taken to keep the system in serviceable 
condition or to bring it back to serviceability. There are two types of maintenance, 
corrective and preventive maintenance [Ref. 5: pp. 1-8]. Corrective maintenance is used 
after a failure. This does not necessarily mean that such action has not been foreseen. 
Preventive maintenance aims to reduce the probability of failure and includes: 

-Planned (scheduled) maintenance in which specified components are replaced 
(usually at regular intervals). The maintenance time is usually based on component 


lifetime distributions. 


-Unplanned (condition-based) maintenance in which the decision to replace or not 
to replace is made according to the outcome of a diagnostic study. 

The policies that are designed to reduce the number or the probability of system 
failure and maintenance costs are called maintenance policies [Ref. 21: pp. 1-3]. We 
concentrate on age replacement policies where a system is replaced at age T or failure, 
whichever comes first. For this problem to make sense, C, must be less than C, and the 
system must age with time i.e. the failure rate of the system must be increasing with age. 
It is not wise to replace the equipment before failure when the system has a constant or 
decreasing failure rate such as when the system failures occur according to an exponential 
distribution. In this case, replacement will not reduce the probability of system failure 
occurring in the next instance of time. When using an age replacement policy the 
question always asked is, “At what age should the equipment be replaced". If the 
replacement occurs too frequently the maintenance costs will be high. If replacement is 
too infrequent, the system will fail more often than necessary and again, the maintenance 
cost will be too high. Thus it is desirable to find an "optimal" replacement age. Here, 
the optimal age is defined as the age which yields the minimum long mun expected 
maintenance costs per unit time. 

We will assume that the maintenance action returns the system to the good-as-new 
condition, thus the same services are provided as before replacement. To accomplish this 
we assume that the systems used for replacement have lifetimes that are independent and 
identically distributed according to a phase-type distribution. The class of phase type 


distributions is large and includes Exponential and Gamma distributions along with 


convolutions and mixtures of these distributions. This fact, along with the fact that phase 
type distributions have a physical interpretation make them particularly well suited for 
modeling system lifetimes. 

The optimal replacement age depends on the underlying phase type distribution. 
This is in general unknown and must be estimated. Estimation for phase type 
distributions based on iid lifetimes has been studied by Neuts and Meier (1980). 
Estimation of the optimal age of replacement for phase type distributions is new. 
However both parametric and nonparametric estimation based on iid observations have 
been considered by Arunkumar (1972), Ingram and Schaeffer (1976), Bergman (1977), 
and Barlow (1978). Here we use a sequential approach for estimating the optimal 
replacement age. At each replacement the estimate is updated and the new system is 
subject to an age replacement policy based on the optimal age estimated so far. This type 
of sequential approach has been studied parametrically by Oclay (1990) and 
nonparametricly by Bather (1977), Frees and Ruppert (1985), Aras and Whitaker (1991), 
Wu (1990) and in a decision theoretic framework by Glazebrook, Bailey and Whitaker 
(1991). 

Phase type distributions are described in Chapter II. The sequential estimation 
procedure to estimate the optimal age of replacement is given in Chapter II, and in 
Chapter IV we present results comparing the sequential estimation procedure assuming 
an underlying phase type distribution with a nonparametric procedure. Conclusion and 


recommendation are given in Chapter V. 


Il. PHASE TYPE DISTRIBUTION 


A. GENERAL 
A phase type distribution is defined as the distribution of the time until absorption 
in a finite state Markov process. To determine the distribution of the absorption time, 


consider an m+1 state, continuous time Markov chain whose infinitesimal generator Q 


has the form 


o- Ea 


where T is a nonsingular m X m matrix with negative diagonal elements and nonnegative 
off-diagonal elements, T° is vector of length m with nonnegative elements and 
0=(0,0,...,0). Moreover, 

Te + T° = 0’ 
where e = (1,1,...,1)’. 

By the construction of Q the state m+1 is absorbing and all other states are 
transient. The necessary and sufficient condition for this is that the inverse T’ exists 
[Ref. 15: pp. 446]. Let the initial probability vector of the Markov chain be given by ( 
Q@, Qn+, ) With we + a,,, = 1 and @ being the m dimensional initial probability vector 
of transient states such that 0 < awe s 1. Let X be the time until absorption into the 
(m+1)* state. The probability distribution F(x) of the time until absorption in the state 


m-+1 corresponding to the initial probability vector (a, a4; ) 18 given by 


F(x) = 1 - aexp(Tx)e, forx 20, (2.1) 
where exp(A) is the matrix exponential of the matrix A defined in Ref. 8. 

The probability distribution F(x) on [0, o) is said to be of phase type (PH- 
distribution) and the pair (a, T) is called a representation of F(x). If a,,, > O the 
distribution F(x) has a jump of height a,,,, at x = 0. 

On (0, ©), the probability mass is distributed according to a density given by 
f(x) = awexp(Tx)T° , for x > 0, (2.2) 
The noncentral moments p; = E[X'] of F(x) are all finite and given by 


u; = (-1)il(eT ie) , for i = 1 [Ref. 15: pp. 446] . (2.3) 


B. COST FUNCTION 

Let X,, X,,... be a sequence of independent and identically distributed (11d) positive 
random variables with distribution function F. The sequence X,, X,,... represents the 
sequence of system lifetimes that would be observable if the system were replaced at 
failure. Let C,, C, (C, > C,) and ¢@ be the respective cost of planned replacement; 
unplanned replacement and the replacement age. The long run expected cost per unit time 


can be verified [Ref. 1: pp. 87] to be 


C.F CF 
i, F(u) du 


where F(u) = 1 - F(u) is the survival function. A sufficient condition that guarantees a 


unique and finite @° that minimizes C(¢), is that F have strictly increasing failure rate. 


When the distribution F is phase type with representation (a, T), the long run expected 


cost function is 


C,[1-aexp (TM) e] + C,[aexp (1) e] 


C(o) 
[Peexp (tu) edu (2 =) 
C, - aexp (T) e[C, -C,] (2.6) 


aT 4[exp(T) - Ije 
where I is an identity matrix. 
The phase type distribution does not necessarily have increasing failure rate, so the 
cost function C(@) does not necessarily have a unique finite minimum. In the m=3 case 
even when the initial probability of the transient states is fixed the minimum can be 


unique and finite or can occur at infinity as shown in Figure 1, where a = (1, 0, 0) and 


-15 13 2 -20 10 10 
fT =|3 -15 12| , T, =|10 -20 10 
2 5 -15 10 10 -25 


In Figure 2, the same phenomena is shown when F has the same nonsingular 3 X3 matrix 


-18 13 5 
T,=)10 -25 15 
4 4 -24 


but the initial probabilities are different. 
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Figure 1 The long run expected cost per unit time curves of PH distribution with 


T, or T, 


the same a@ and different T 
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LONG RUN EXPECTED COST PER UNIT TIME OF PH DISIN. 
WITH REPRESENTATION T3 AND C1=1, C2=5 


eee ALPHA = (0.5, 0.3, 0.2) 
ALPHA = (1.0, 0.0, 0.0) 





1 2 2 
REPLACEMENT AGE 


Figure 2 The long run expected cost per unit time curves of PH distribution with 


the same T and different a 


I. THE SEQUENTIAL ESTIMATION PROCEDURE 

In theory, any distribution of a nonnegative random variable can be approximated 
arbitrarily closely by a PH distribution [Ref. 11: pp. 1-2]. This means that the PH 
distributions are dense in ¥ where F¥ is the set of distributions with support on [0, °°). 
The paper of Johnson and Taaffe (1988) shows that the finite mixtures of Erlang 
distributions and PH distributions with a Coxian representation are both dense in ¥. Since 
both families are subsets of the family of PH distributions, this implies that PH 
distribution are dense in ¥ [Ref. 10: pp. 1-8]. Thus, the PH distribution can be used to 
approximate any unknown distribution of lifetimes (lifetime is always greater than or 
equal 0). Another feature that makes PH distributions a desirable choice to model system 
lifetimes is that they may have a physical interpretation. The absorbing state represents 
system failure and the transient states may represent different levels of a system’s ability 


to function. 


A. PARAMETER ESTIMATION 

There are many ways to estimate the parameters of a distribution including the 
method of maximum likelihood (MLE), and the method of moments. The PH 
distributions have special properties which make estimation difficult. First, the number 
of parameters is not fixed, it varies with the number of transient states, e.g. a PH 
distribution with 2 transient states has 6 parameters, a PH distribution with 3 transient 


states has 12 parameters etc. Also the probability distribution and density function include 


the exponential of a matrix. This makes it hard to work with the likelihood function. In 
this thesis, the moment matching method is used to estimate the parameters of a phase 
type distribution but the MLE method and its associated problems are also discussed . 


In addition, the number of transient states, m, is assumed to be known. 


1. The Maximum Likelihood Method 
Let X,, X,,... represent the sequence of system lifetimes, where X,, X,,... 
are iid PH distributions with representation (a, T). After N observations the data 
available for estimation will be (Z, 5,),(Z,,6,),...,(Zy,dn) where Z;=min (X,,6°;,), X; is 
the i“ lifetime and $°,, is the optimal age replacement estimated after i-1 observations and 
6,=1 if X,<$',, and 6,=0 otherwise. Assuming system lifetimes have a PH distribution 
the likelihood for the replacement ages Z,,....Zn and types of replacement 6,, 63,...,dy 


iS 


1-8, 


L(e,T) = [aexp (TZ,) T*] *t [aexp (TZ,) e] 


N 
4{=1 


N 
= Il [-awexp (TZ,) Te] *t [wexp (1Z,) e] 
1-1 


Tet 


It is too difficult to maximize this likelihood directly by differentiating L(a,T) and 
solving for the MLE’s. There are several reasons for this. The likelihood and its 
derivatives include the exponential of a matrix form that cannot be written in closed form 
and the parameters are subject to numerous constraints. In addition, there are many 


likelihood equations due to the number of parameters (3 transient states will have up to 
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12 equations). An alternate approach is to find approximate MLE’s using nonlinear 


programming algorithms 1.e.: 


N 
MAX L(a,T) = [[ [-aexp(TZ;) Te] *t [aexp (TZ,) e] 1-84 
4=1 


S.T. a > 0 
ae = 1 
t; < O for1 = 1,2,...,m 
t 2 0 fori ¥ j 
eee = Gs fori = 1,2,...,m 
det(T) + 0 


where m = the number of transient states of the PH distribution. Even when m 1s 
assumed to be known, there is still the problem of approximating the exponential of a 
matrix and the optimization software to take care of this problem is not available. Thus, 


it is very difficult to use this approach to estimate the PH distribution parameters. 


2. |The Moment Matching Method. 

The PH distribution is a complicated distribution, there are many parameters 
and there are difficulties with computing the exponential of a matrix. One property, the 
existence of T’, implies that all noncentral moments of a PH distribution are finite. The 
k* moment is given by 

wm, = (-l¥k!aT’e. (3.1) 
Moment matching is a common method for estimating. Using moment matching by- 


passes the problem of evaluating the exponential of a matrix [Ref. 12: pp.3]. 
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Let m be the number of transient states of the PH distribution. There are m(m+1) 


parameters that need to be estimated, thus the m(m+1) equations from the moments that 


need to be solved are 


By — -aT'e 
in = 2!aT'e 
Finmet) = (-1)™™*?(n’+m)!aT’e, 


where ji, is the k™ sample moment. 

This method is still inappropriate due to the large number of equations and 
the fact that when the dimension of matrix T is big and the moments are of high order, 
substantial amount of error is introduced when trying to solve these equations. Johnson 
and Taaffe have shown that the moment matching method and nonlinear programming 
approaches can be used to estimate the parameters of a PH distribution. They match only 
the second and third standardized moments [Ref. 11: pp. 3-11]. 

Let py; be the i* central moment. The second standardized moment, the coefficient 


of variation (C), defined as 


Coz standard deviation 
mean 


can be written as 
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(p,) a7 2 


C (3.2) 
asl 
The third standardized moment, the coefficient of skewness, (vy) is defined as 
Ht; 
Yo. > sews * (3.34) 
( i, ) Biz 


And the t® standardized moment is p,/(,)”? fort = 3,4,... . 
The reason for matching standardized moments is that the standardized moments do not 
change with scale changes. So, the shape of a PH distribution with representation (a,T) 
will not change if it is multiplied by a nonnegative number. 

The nonlinear programming formulation for matching C and y to a 
continuous PH distribution is given by 


MIN. w,(C-C(0))? + wy. (y-y(0)) 


S.T. 
a = 0 
ae = | 
t; < 0 fori = 1,2,...,m 
t; = 0 fori ¥ j 
Ving ty S -t; fori’ —s24m 
det(T) # 0 
where 


C(O) = the second standardized moment of the current solution. 


(0) = the third standardized moment of the current solution. 
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m = the number of transient states of the PH distribution. 

t; = the elements of row i and column j of nonsingular matrix T. 

W,,W, are positive weights chosen to guide the search and w, = 3w,. 

The moments are matched when the objective function value is zero. The solutions which 
match the target moments are called "moment-matching solutions" [Ref. 11: pp. 5]. 

Even with fixed dimension, this procedure may have many solutions and these 
solutions are unpredictable. There are some practical points that need to be made 
associated with the properties of PH distribution and the problem of using nonlinear 
programming to get the moment matching solutions: 

- Different combinations of initial solutions and algorithm parameters may lead to 
different moment-matching solutions. Also, some combinations of initial solution and 
algorithm parameters may not lead to a moment-matching solution and the nonlinear 
programming algorithm may not converge. 

- For a given dimension we do not know how many solutions exist or how to guide 
the search toward a preferable solution. 

- The moment-matching solutions do not guarantee that the estimated cost function 
will have a finite minimum. When choosing among several solutions, a solution which 
gave a cost function with a finite minimum was always chosen. 

- The values of w, and w, can be modified to guide the search; generally w, =1, 
w,21 and w, = 3w,. The magnitude of w, and w, can be increased to get more or 


sufficient accuracy [Ref. 11: pp. 6]. 
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- If the feasible solutions are known, parameters can be bounded to guide the 
search. However, some bounds on the parameters may not lead to a moment-matching 
solution or to a solution whose cost function does not have a unique finite minimum. 

- User interaction is often necessary in obtaining the appropriate or preferable 


solutions. 


B. THE SEQUENTIAL ESTIMATION PROCEDURE 

Let X,,X,,.....,Xy be the sequence of lifetimes of the system and {#',} be the 
sequence of estimators of ¢° where $°, is the estimator after N replacements. Then after 
N observations, we have the right censored data (Z,,6,),(Z,,6),....,;(Zy,6n)- 
where Z, = min (X%,,$';.,) 

6 =1 if X; s $;,, otherwise 6, =0. 

Because the data are right censored the usual method of moments approach for estimation 
needs to be modified. Rather than use sample moments calculated from an empirical 
distribution, we use moments from a nonparametric estimator of F based on the censored 
data. From the right censored data after N observations, the procedure to compute the 
estimators {°;} is developed as follows. 

1. Use the right censored data sequence (Z,,6,),(Z),6,), ...,(Zy,6x) to estimate 


moments with F=1-F by the product-limit estimator: 
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F(t) = toe (3.4) 
(112.4) <0) N-i+1 





where Z,), Z,.--, Zqy are the order statistics of Z,, Z,,..., Zy and by), @),---, San are 
ordered according to the ordering of Z,), Zg,..., Zqy. Then calculate probabilities P,;, 
associated with Z,, by 

P, = F(Z) : F(Z) (3.5) 
Note, since F has discontinuities only at Z~ where 6, = 1, that P, = 0 when 6, = 0. 


Then estimates of the k™ moment are given by 
~ N 

EB(X*) = YO PZ 3) (3.6) 
ie1 


The estimate of the 2"? standardized moment is 


VE(X*) - (Elda (3.7) 
E(X] 


> 


and the estimate of the 3" standardized moment is 


E(X?) -3B[X] BE[X?] +2 (B(x) )? 


(VE[X?] - (B[X] )?)3 4 


2. Estimate parameters of PH distribution (a@,T), by matching the second and third 


(3.8) 


standardized moments using the nonlinear programming approach. 
When executing the nonlinear programming algorithm using the model described 
in the previous section, the det(T) + 0 constraint is replaced by a constraint on the 


expected lifetime. The reason for doing this is that we could not formulate an algebraic 
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constraint equivalent to the constraint det(T) ~ 0 . The problem can be solved by a 


constraint on the expected lifetime, i.e. by taking 
afte = P[X] (3.9) 


where E[X] is calculated in step 1. This constraint will make the search for preferable 
solutions easier. The initial solution is chosen by the user. Some initial solutions may or 
may not lead to a moment matching solution. Typically, the initial vector a@ is 
(1/m,1/m,...,1/m) or (1,0,...,0) and an initial matrix T consists of t;; = -1 andt,; = 1/m 
fori * j, i = 1,2,...,m. In this thesis the initial vector @ is (1,0,...,0) and an initial 
matrix T is taken to have t;; = -1 and t,; = 1/m (Ref. 11: pp. 9] . 

3. Using the parameters estimated in step 2, the estimated cost function is 


C, - &exp (TM) e[c, -C] 


C(d) = = 
aT * [exp (I) - Ile 


(3.10) 


The new estimate of the optimal age of replacement $°, is taken to be the ¢ that 
minimizes (3.10) by enumerative search. 

4. Compare the $°, with X,,, then repeat Step 1. 
The initial solutions for matrix T and vector «& in step 2 are taken to be i previously 
estimated values of T and a. In case the moment matching solution cannot be met, the 
initial solution in step 2 will be taken to be the original initial solution and the previous 
optimal age replacement is used for step 3. 


The replacement cost for i* system is C, if X; < $°;,; otherwise the replacement 


cost is C,. The actual total replacement cost for the first N systems that are observed is 
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N 
Cy = py, | Grr ee erachaGrs (3.11) 


1=1 


and the total operating time for the N systems is 
N 
Cy 7" > Min gir $* ,_.) = >> 2; e (3.12) 
i=1 
So, the actual average cost (AAC) after N replacement is computed by 
AAC, = — (3.13) 


In this thesis, the moment matching nonlinear programming approach for estimating 
parameters of PH distributions uses the GAMS program as shown in Appendix A (for 
a PH distribution with 3 transient states). The Fortran programs are written to estimate 
the optimal age of replacement, second and third standardized moments and for data 


preparation as shown in Appendix B. 
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IV.COMPARISON WITH NONPARAMETRIC PROCEDURE 


A. NONPARAMETRIC PROCEDURE 

An alternate approach for estimating ¢° is to estimate nonparametrically, as in Aras 
and Whitaker (1991). The procedure is much simpler computationally than the parametric 
procedure described in the previous section for PH distributions. In the nonparametric 
procedure after N replacements the product limit estimator F of F given in (3.4) is used 


to estimate the cost function C(¢) as 


e C.P(b) + CF(O) 
Cylo) = a lil (4.1) 
f F(u) du 


where F = 1 - E is the estimator of the cdf F, the estimator °, of ¢ is then found by 
minimizing C,(¢). In general one would expect parametric procedures to do much better 
than nonparametric procedures. However, with the numerical difficulties involved with 
estimating parameters for PH distribution and the fact that the family of PH distributions 
is so large it is not obvious which procedure is best. The criterion for Porpanicen is the 


actual average cost per unit time, AAC,, after N replacements. 
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B. COMPARISON OF NONPARAMETRIC AND PARAMETRIC PROCEDURE 
FOR THE PH DISTRIBUTION 
Simulation was used to compare sequential estimation based on PH distributions 
with nonparametric sequential estimation, C, = 100, C, = 500. System lifetimes were 
generated by an increasing failure rate PH distribution with representation (T,, ,a,;) 


where, 


-5 4 1 
T,=|1 -6 5| , @, = (10090 on cme) 
1 1-7 


For this representation the long run expected cost per unit time is given in Figure 3. 

To give the parametric procedure a chance, the first 25 system lifetimes are 
uncensored. After the first 25 observations, the parameters of the PH distribution are 
estimated sequentially as in Chapter III Section B. The actual average costs (AAC) are 
calculated and compared for small, medium and large sample sizes with planned cost C, 
= 100 and unplanned cost C, = 500. The 20 initial data sets are simulated as shown in 
the Appendix A. 

The programs are separated into 5 parts. The first one is written in GAMS as 
shown in Appendix A, to estimate the parameters of a PH distribution by the nonlinear 
programming approach. The rest are written in Fortran as shown in Appendix B, to 
estimate the optimal age of replacement based on the estimated PH distribution from the 


GAMS program; to simulate the system lifetime; to prepare the right censored data; to 
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estimate the expected lifetime, second and third standardized moments, and finally to 


tabulate the results and to prepare the initial data for the next estimation. 


LONG RUN EXPECTED COST PER UNIT TIME OF PH DISTN. 
WITH REPRESENTATION 141, a41 AND C1=100, C2=500 
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Figure 3 Long mun expected cost per unit time curve of PH distribution with 
representation T,, and a, 


In this simulation, the number of transient states is fixed at 3 and the initial 
probability vector a = (1, 0, 0). By fixing a = (1, 0, 0) we reduce the number of 
parameters to be estimated and use a model in which all systems start in the same state, 
when states are thought of as the level of system repair, this choice better reflects the 


idea that replacement systems are as good as new. Bounded parameters and user 


21 


interaction are used in guiding each estimation during the simulation to a preferable 
solution. 

Tables 1, 2 and 3 summarize the simulation results of the actual average cost 
resulting from sequential estimation of PH and the nonparametric procedure. Also given 
are the signed ranks of the differences in the actual average costs. These are used in the 
Wilcoxon Signed-Rank test to test the hypothesis: 

H,: E[AACpy] = E[AACyonp] 

H,: E[AAC,,] < E[AACyonp] 
by using T* the sum of the positive ranks as the test statistic . 
Tables 1, 2 and 3 give the statistic T* as 35, 10 and 2 respectively. When compared to 
the value in Table 9 [Ref. 13: pp. 780] at N=20 all p-values are less than 0.005 which 
implies that the sequential estimation of a PH distribution is better than sequential 
estimation of nonparametric for all sample sizes (small, medium and large). The box 
plots in Figures 5 and 6 of actual average cost versus the number of replacements at 
different sample sizes show that the actual average costs decrease as the number of 
replacements increase and that AAC,,, for the parametric procedure decrease more than 
the AACyonp for the nonparametric procedure. 

A second, PH distribution was chosen to compare the parametric PH procedure 
with the nonparametric procedure. The second PH distribution has an average long run 
expected cost function that is more shallow than the first PH distribution. It has 


representation T,, and a,, 


is 


where 


5 4 1 
rr. eos, Clty, = «(1.0,0.0,0.0) , 
1 2 -7 


and the long run expected cost per unit time is given in Figure 4. 


LONG RUN EXPECTED COST PER UNIT TIME OF PH DISTN. 
WITH REPRESENTATION 142, a42 AND Ci=100, C2=500 
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Figure 4 Long run expected cost per unit time curve of PH distribution with 
representation T,, and a,, 


Tables 4, 5 and 6 summarize the result and give the statistic T* of 93, 81 and 47 


respectively. The p-value taken from table 9 [Ref. 13: pp. 780] indicated that the 
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parametric procedure is not better than the nonparametric procedure for small and 
moderate sample sizes. In fact, for some cases (e.g. runs 16 and 20 in Table 4 and rns 
1, 5 and 16 in Table 5) the actual average cost resulting from using the nonparametric 
procedure can be quite a bit less than the actual average cost resulting from the 
parametric procedure. This is due to the fact that C(¢) is shallow at ¢°, thus even though 
C(¢) may be a reasonable estimator of C(¢), the variance of its minimizing value $° is 
higher than for the previous PH distribution. This means that ¢° is often underestimated. 
Because C(¢) increases so rapidly to the left of ¢°, underestimating ¢° can increase the 
actual average costs dramatically. Another difficulty is that the estimated C(¢) may be 
relatively flat around ¢°. This causes numerical difficulties with the nonlinear 
programming algorithm. For large sample sizes the parametric procedure does do better 
than the nonparametric procedure. Here C(@) estimates C(#) more closely and it is easier 
for the nonlinear programming algorithm to identify the correct solution. 

In Figures 5 and 6, the improvement in actual average cost with sample size is 
shown for both the parametric and nonparametric procedure for the PH distribution with 
representation T,, and a,,. As expected, both procedures improve (actual average costs 
decrease) with sample size. However, the actual average cost for the parametric 
procedure decreases faster than for the nonparametric procedure. Both procedures exhibit 
considerable variability. The solutions have got still mixed that we should do more study 


on another representation until the more precise solutions come up. 
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Table 1 COMPARISON OF ACTUAL AVERAGE COSTS FOR SMALL SAMPLE 
SIZES (N = 50) OF PH DISTRIBUTION WITH T,, AND og, 


| 1 [90.095 | 754.555 | 759.142 | 64.45 |e 
5_| 74330 | 757.528 | 735.0588 | 43.198 | 8 
7 
10 
| sia99 | 767.586 | so1.tso | 47.403 | 9 
2 | 621 | 667.553 |or.16 | 22.632 | 
13 
e29.204 | 667.805 | 682.719 | 38st | 6 

688.572 
6_ | 636.963 


583.226 625.161 609.036$ | -41.935 7 


644.065 654.691 676.005 -10.626 


AAC ,, = Actual average cost replacing when failure occurs 
DIFFERENCE = AAC), - AACyonp 
* positive rank 
Tt = 35; 
p-value < 0.005 
$ AACg, is less than only AACyonp 
! AAC, is less than AAC,,, and AACyonp 


— — 
~~) 
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Table 2 COMPARISON OF ACTUAL AVERAGE COSTS FOR MEDIUM 
SAMPLE SIZES (N = 100) OF PH DISTRIBUTION WITH T,, AND a,, 


Run || AACm | AAC | _AACye_| DIFFERENCE | RANK _ 
£90,971 
649.886 
96.858 
90.195 


77.742 607.414 610.850 -29.672 





bo 


ms 
ON 


’ 2 
: 5 

86.414 766.365 762.245 $ “19.951 

50.106 612.046 638.286 -61.94 


~] 


a ee ee ee ee ee eee ee : 
JT Tnlmai Pp [wo frm |= [ao ® Wo 


7 


37.268 696.249 689.552 $ -58.981 
84.253 730.138 786.920 -45.885 


a - NO 
% ¥ ¥ 


02.039 
15.087 
77.302 
86.838 
26.213 
26.605 
71.908 
32.954 


oo 


0 
11 
17 
12 
16 
17 
3 
13 
) 


9 || 635.989 679.399 | 667.1419 | -43.41 9 
20 || 590.881 657.064 701.735 66.183 15 


AACpgr = Actual average cost replacing when failure occurs 
DIFFERENCE = AAC), - AACyonp 
* positive rank 
T* = 10 ; 
p-value < 0.005 
$ AACgz is less than only AACyonp 
! AACgz is less than AAC,, and AACyonp 
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Table 3 COMPARISON OF ACTUAL AVERAGE COSTS FOR LARGE SAMPLE 
SIZES (N = 200) OF PH DISTRIBUTION WITH T,, AND a, 





cd ee ee ee 


Li | 622.560 | 6ss.eo1 | oss.s6ss | -33.301 | 7 
| 2 | 6se.cso_ | 11.831 | roasoss | ssn | on 
4 | 67.141 | 749.866 | 754.090 | mms |e 
| 5 | erases | r14201 | ro7srs | 39.03 | 9 
| 8 | sas.ess | 61s.s20 | 647.927 | 604s | 13 
| 8 | sam | 607.068 | 649.386 | 25.357 |S 
10 _| 655.709 | 717.569 | 766.000 | -oi.86 | a 
12 | sao.oss__| 649.545 | 81.004 | -60.457_ | 18 
13 | 602.032 | 643.703 | ao.ss2_ | 478 | 10 
-is_| 685.977 | 767.133 | 755.3318 | tse | 18 
-1s_| 646.752 | 3.719 | 663.498 | 76.967 |e 
17 | 80.457 | 713.146 | 754.616 | -32.689 | 6 
SS ee 

| 636.649 | 74.845 | o7o.ises | 38.196 | 8 


636.649 674.845 670.146 $ -38.196 


588.468 | 654.994 670.194 66526 | 15 


AAC ,, = Actual average cost replacing when failure occurs 
DIFFERENCE = AAC,,, - AACyonp 
* positive rank 
T* =2 ; 
p-value < 0.005 
$ AACp, is less than only AACyonp 
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Table 4 COMPARISON OF ACTUAL AVERAGE COSTS FOR SMALL SAMPLE 
SIZES (N = 50) OF PH DISTRIBUTION WITH T,, AND ay 


098 -80.466 19 
no il teal 47 
mim ile 
rm Per falas 
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ej po 
29.150 -39.118 
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577.580 576.101 600.152 1.479 


560.629 519.418 532.841 41.211 


AACp, = Actual average cost replacing when failure occurs 
DIFFERENCE = AAC,,, - AACyonp 
* positive rank 
Tt = 93; 
p-value > 0.05 
$ AACgp, is less than only AACyonp 
! AACgs iS less than AAC,, and AACyonp 
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Table 5 COMPARISON OF ACTUAL AVERAGE COSTS FOR MEDIUM 
SAMPLE SIZES (N = 100) OF PH DISTRIBUTION WITH T,, AND a, 


bo 


Le) —_ — poe rad fom — Pees a po ret 


75.157 
02.762 
35.426 
80.480 


19 * 
4 
20 
Z 
7 
1 


25.770 
18.669 
34.978 
75.021 18 * 
35.911 17 
17.194 | 636.912 | was.7s | 19.718 


oT Fa aE Snr. hr 
539.256 | 520.034 531.189 19.222 ise | 


AAC,, = Actual average cost replacing when failure occurs 
DIFFERENCE = AAC,,, - AACyonp 
* positive rank 
Tt = 81 ; 
p-value > 0.05 
$ AACy, is less than only AACyonp 
! AAC), is less than AAC,,, and AACyonp 
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Table 6 COMPARISON OF ACTUAL AVERAGE COSTS FOR LARGE SAMPLE 
SIZES (N = 200) OF PH DISTRIBUTION WITH T,, AND ay 


578.172 










594.573 579.915 588.547 14.658 
580.261 551.976 560.463 28.285 
526.351 558.604 573.294 -32.253 


an 
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632.238 624.275 
612.338 
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52.09 


578.172 577.936 607.344 0.236 1 


ok 
20 538.683 530.002 570.516 8.681 4 * 


AACzr = Actual average cost replacing when failure occurs 
DIFFERENCE = AACp, - AACyonp 
* positive rank 
T* = 47 ; 
0.01 < p-value < 0.025 
$ AACg, is less than only AACyonp 
! AAC), is less than AAC,,, and AACyonp 
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Figure 5 The box plot of actual average cost for different sample sizes of 
phase type sequential procedure for T,, and a4; 
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Figure 6 The box plot of actual average cost for different sample sizes of 
nonparametric sequential procedure for T,, and a, 
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V.CONCLUSIONS AND RECOMMENDATIONS 

In this thesis, the age replacement policy has been considered. A system is replaced 
at the time of failure or at a scheduled replacement age whichever comes first. The 
replacement system is assumed to be as good as new. The objective is to achieve the 
minimum long run expected maintenance cost per unit time. The system lifetimes are 
assumed to have PH distributions. Estimating the parameters of a PH distribution is 
difficult because it may have many parameters and its density function involves a matrix 
exponential. Moment matching with a nonlinear programming approach is used for 
estimating the parameters. This method does not guarantee that the estimated cost 
function will have a unique and finite minimum. In this thesis we restric our attention to 
the case with initial probability vector @ = (1, 0, 0). Sequential estimation using the 
phase type procedure gives smaller average costs than the nonparametric procedure for 
both small and large sample sizes when the underlying long run average cost function has 
a very distinct minimum at ¢°. When this cost function is shallow around ¢’, the 
parametric procedure does not do better than the nonparametric procedure for small 
samples. In fact, for cost functions that are shallow, it is better not to use a maintenance 
policy i.e. to take $° = o, until the sample sizes are large enough to give reasonable 
estimates of the underlying parameters. The reason for this is that for such cost 
functions, overestimating ¢° does not increase the long run average cost much, at the 


same time the decrease in the amount of censoring with over estimating ¢° greatly 
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improves estimation of the underlying distribution function F. Conversely 
underestimating ¢° for these cost functions causes a drastic increase in long run average 
costs and increase censoring making estimation very difficult. 

Recommendation for the future research 
The following forms a list of future research initiatives: 

-Examine the impact of the product-limit estimator in estimating standardized 
moments on the rest of procedure. 

-Do more experimentation on different representations of phase type distributions 
to get more precise conclusion. 

-Look for other policies using phase-sensitive estimate and modified replacement 
costs over time. 

-The phase-type sequential estimations we have been so far assume the underlying 
lifetimes are iid and after replacement the system is new. To be more realistic, we can 
modify the initial probability vector, increase transition rates between states, or both, 
over time. 

-Use phase-type sequential estimation when the underlying life distribution F comes 
from Gamma, Weibull etc. that have increasing failure rate and compare with the 
nonparametric procedure. Since we can use the denseness of phase type distribution 
property to approximate the set of distribution with support on [0, a) 1.e., the set of 


lifetimes, we may get a better method when the distribution of lifetimes is not known. 
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APPENDIX A. 


$TITLE PH DISTRIBUTION ESTIMATION PARAMETERS 
a 


*_---------- GAMS AND DOLLAR CONTROL OPTIONS -------------------------------------- 
- (See Appendices B & C) 
OPTIONS 


WORK = 100000, 
LIMCOL = 0 , LIMROW = 0, SOLPRINT = OFF , DECIMALS = 2 
RESLIM = 100, ITERLIM =100000, OPTCR = 0.0 , SEED = 3141; 


* This program uses nonlinear programming to estimate the 3-transient states of phase 
* type distribution parameters. Matching the second and third standardized moments are 
* used. All adjustable data are prepared by the "FILE SCALAR2" and uses the 

* "SINCLUDE" statement bring to the program. 


$INCLUDE "FILE SCALAR2" 


VARIABLE 

objective function value 

first entry of the matrix 

second entry of the matrix 

third entry of the matrix 

fourth entry of the matrix 

fifth entry of the matrix 

sixth entry of the matrix 

seventh entry of the matrix 

eighth entry of the matrix 

ninth entry of the matrix 

ALI initial probability of being in state 1 
AL2 initial probability of being in state 2 
AL3 initial probability of being in state 3 ; 


srmantMOaAWPSN 


EQUATION 
OBJ Define objective function 
EQ1 Set initial probability of being in state 1(nonnegative) 
EQ2 Set initial probability of being in state 2(nonnegative) 
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EQ3 Set initial probability of being in state 3(nonnegative) 

EQ4 Sum over all of initial prob. less than or equal 1.0 

EQS Set the first diagonal entry to be nonpositive 

EQ6 Set the second diagonal entry to be nonpositive 

EQ7 Set the third diagonal entry to be nonpositive 

EQ8 Set off diagonal entries in the first row to be positive 

EQ9 Set off diagonal entries in the first row to be positive 

EQ10 Set sum off diagonal of the first row less than or equal first diagonal 
EQ11 Set off diagonal entries in the second row to be positive 

EQ12 Set off diagonal entries in the second row to be positive 

EQ13 Set sum off diagonal of the second row less than or equal second diagonal 
EQ14 Set off diagonal entries in the third row to be positive 

EQ15 Set off diagonal entries in the third row to be positive 

EQ16 Set sum off diagonal of the third row less than or equal third diagonal 
EQ17 Define the expected lifetime ; 


* MINIMIZE 
OBJ.. 
7 =E= 


#«* 3%SQR(C - C(0)) *** 


3*POWER((((2*AL1 “(POWER ((-(F*H) + E*D),2) +(C*H-B*I)*(F*G-D*I) 

+(-(C*E) +B*F)*(-(E*G) + D*H) 
+(-(F*H) + E*1)*(C*H-B*D + (C*H-B*I)*(-(C*G) +A*I) + (-(C*E) +B*F)*(B*G-A*H) 
+(-(F*H) + B*1)*(-(C*E)+B*F) +(C*H-B*I) *(C*D-A*F) + (-(C*E) +B*F)*(-(B*D) + 
A*B)) 


+ 2*AL2*((F*G-D*I)*(-(F*H) +E*D) +(-(C*G)+A*l)*(F*G-D*I) +(-(E*G)+D*H)* 
(C*D-A*F) 

+(F*G-D*I)*(C*H-B*I) + POWER((-(C*G)+A*D),2)+(C*D-A*F)*(B*G-A*H) 

+(F*G-D*I)*(-(C*E) +B*F) + (-(C*G) +A*I)*(C*D-A*F) +(C*D-A*F)*(-(B*D)+A*B)) 


+ 2*AL3*((-(B*G)+D*H)*(-(F*H) +E*l +(B*G-A*H)*(F*G-D*I + (-(B*D)+A*E)* 

(-(E*G) + D*H) 
+(-(B*G)+D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*D+(-(B*D) +A*E)* 
(B*G-A*H) + (-(E*G) + D*H)*(-(C*E) +B*F) + (B*G-A*H)*(C*D-A*F) + 
POWER((-(B*D)+A*B),2)))/ 


POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*D,2) 
-POWER((-(AL1*(-F*H+E*I+C*H-B*I-(C*B) +B*F) 
+AL2*(B*G-D“I-C*G+ A*I+C*D-A*F)+AL3*(-B*G+D*H+B*G-A*H-B*D+A*B))/ 
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(-(C*E*G) + B*F*G + C*D*H-A*F*H-B*D*I+ A*E*D),2))**0.5 / 
(-(AL1 *(-(F*H) + B*I+C*H-B*I-(C*E)+B*F) 

+ AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 
+AL3*(-E*G+D*H+B*G-A*H-B*D+A*B))/ 

(-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I)) - MM2),2) 
*** W2*SQR(GAMMA - GAMMA(0)) *** 


+ POWER(((((-6*AL1*((POWER((-(F*H) +E*I),2)+(C*H-B*D *(F*G-D*1) + 
(-(C*E) + B*F)*(-(E*G) + D*H))*(-(F*H) +E*T) 


+(((F*H) + E*1)*(C*H-B"1) + (C*H-B*I)*(-(C*G) +A") + (-(C*E)+B*F)* (B*G-A*H))* 
(F*G-D*I) 


+((-(F*H) + E*I)*(-(C*E) + B*F) + (C*H-B*I)*(C*D-A*F) + (-(C*E) + B*F)* 
(-(B™D) + A*E))*(-(B*G) + D*H) 


—-+—=—-—~S~«~—iS;CS:CS;7SC SF +B*F)* (-(E*G)+D*H))* 
(C*H-B* 


+((-(F"H) + B*1)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I) + (-(C*E) + B*F)*(B*G-A*H))* 
(-(C*G)+A*I) 


+ ((-(F*H) +E*I)*(-(C*E)+B*F)+ (C*H-B*I)*(C*D-A*F) + (-(C*E)+B*F)* 
(-(B*D)+A*E))*(B*G-A*H) 


+ (POWER ((-(F*H)+E*I),2)+(C*H-B*I)*(F*G-D*1) + (-(C*E)+B*F)* 
(-E*G) +D*H))*(-(C*E) +B*F) 


+ ((-(F*H) + E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*T) + (-(C*E)+B*F)* 
(B*G-A*H))*(C*D-A*F) 


+ ((-(F*H) +E*D)*(-(C*E)+B*F) + (C*H-B*I)*(C*D-A*F) + (-(C*E)+B*F)* 
(-(B*D) + A*E))*(-(B*D) + A*E)) 


-6*AL2*(((F*G-D*1)*(-(F*H) + BD) + (-(C*G)+A*1)*(F*G-D*1) + (-(B*G)+D*H)* 
(C*D-A*F))*(-(F*H) +E*D) 


+((F*G-D*I)*(C*H-B*I) + POWER((-(C*G)+A*D),2)+(C*D-A*F)*(B*G-A*H))* 
(F*G-D*I) 


+((F*G-D*I)*(-(C*E) +B*F) + (-(C*G)+A*I*(C*D-A*F) + (C*D-A*F)* 
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(-(B*D) + A*E))*(-(E*G) + D*H) 


+ ((F*G-D"T)*(-(F*H) + E*D) + (-(C*G) + A*T)*(F*G-D*I) + (-(E*G) + D*H)* 
(C*D-A*F))*(C*H-B*D 


+((E*G-D*I)*(C*H-B*I) + POWER((-(C*G) + A*I),2)+(C*D-A*F)*(B*G-A*H))* 
(-(C*G) + A*T) 


+((E*G-D*I)*(-(C*E) +B*F)+ (-(C*G)+ A*I)*(C*D-A*F) +(C*D-A*F)* 
(-(B*D) + A*E))*(B*G-A*H) 


+ ((F*G-D*1)*(-(F*H) + E*D +(-(C*G) + A*1)*(F*G-D*I) + (-(E*G) + D*H)* 
(C*D-A*F))*(-(C*E) + B*F) 


+((F*G-D*I)*(C*H-B*I) + POWER((-(C*G) + A*I),2)+(C*D-A*F)*(B*G-A*H))* 
(C*D-A*F) 


+((E*G-D*I)*(-(C*E) + B*F) + (-(C*G)+A*D*(C*D-A*F) + (C*D-A*F)* 
(-(B*D) + A*E))*(-(B*D)+A*E)) 


-6*AL3*(((-(E*G) + D*H)*(-(F*H) + E*l) + (B*G-A*H)*(F*G-D*I) + (-(B*D) + A*E)* 
(-(E*G) + D*H)*(-(F*H) +E*D) 


+((-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*1) + (-(B*D) + A*E)* 
(B*G-A*H))*(F*G-D*"T) 


+ ((-(E*G) -+D*H)*(-(C*E)+B*F) + (B*G-A*H)*(C*D-A*F)+ 
POWER((-(B*D)+A*B),2))*(-(E*G) +D*H) 


+ ((-(E*G) + D*H)*(-(F*H) + E*D) + (B*G-A*H)*(F*G-D*T) + ((B*D) + A*E)* 
(-(E*G) + D*H))*(C*H-B*I) 


+ ((-(E*G) + D*H)*(C*H-B*T) + (B*G-A*H)*(-(C*G)+A*D) + (-(B*D)+A*E)* 
(B*G-A*H))*(-(C*G) + A*D) 


+ ((-(E*G) +D*H)*(-(C*E)+B*F) + (B*G-A*H)*(C*D-A*F) + 
POWER ((-(B*D) + A*E),2))*(B*G-A*H) 


+((-(B*G) + D*H)*(-(F*H) + ED) + (B*G-A*H)*(F*G-D*I) + (-(B*D) +A*E)* 
(-(E*G) + D*H))*(-(C*E) +B*F) 


+ ((-(B*G) + D*H)*(C*H-B*1) + (B*G-A*H)*(-(C*G) + A*T) +(-(B*D)+ A*B)* 
(B*G-A*H))*(C*D-A*F) 
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+((-(E*G) + D*H)*(-(C*E) + B*F) + (B*G-A*H)*(C*D-A*F) + 
POWER((-(B*D) + A*E),2))*(-(B*D) + A*E))) 


/POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I),3) 


+3*((AL1 *(-(F*H)+EB*I+C*H-B*I-(C*E)+B*F)+AL2*(F*G-D*I-C*G+A*I+C*D- 
A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*B))/(-C*E*G+B*F*G+C*D*H-A*F 
*H-B*D*I+A*B*]))* 


*** expected of lifetime square *** 


((2*AL1 *(POWER((-(F*H) + E*I),2)+(C*H-B*I) *(E*G-D*I) + (-(C*E) +B*F)* 
(-(E*G) + D*H) + (-(F*H) + E*I)*(C*H-B*I) +(C*H-B*I)*(-(C*G) + A*D) + 

(-(C*E) +B*F)*(B*G-A*H) + (-(F*H) + E*I)*(-(C*E)+B*F) +(C*H-B*I) *(C*D-A*F) 
+(-(C*E)+B*F)*(-(B*D)+A*B)) 


+2*AL2*((F*G-D*I)*(-(E*H) +E*I) +(-(C*G)+A*I)*(E*G-D*1) 
+(-(E*G) +D*H)*(C*D-A*F) +(F*G-D*I)*(C*H-B*I) +POWER((-(C*G) +A*D),2) 
+(C*D-A*F)*(B*G-A*H) + (F*G-D*I)*(-(C*E) +B*F) + (-(C*G)+A*I)*(C*D-A*F) 
+(C*D-A*F)*(-(B*D)+A*B)) 


+2*AL3*((-(E*G) +D*H)*(-(F*H) +E*l) + (B*G-A*H)*(F*G-D*1) 
+(-(B*D)+A*E)*(-(E*G)+D*H) + (-(E*G) + D*H)*(C*H-B*1) + (B*G-A*H)* 
(-(C*G) + A*1) +(-(B*D) + A*E)*(B*G-A*H) + (-(E*G) + D*H) *(-(C*E) +B*F) 
+(B*G-A*H)*(C*D-A*F) + POWER((-(B*D) + A*E),2)))/ 


POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I),2)) 
*** 2 time cube of the expected lifetime *** 


-2*POWER(((ALI *(-(F*H) + E*I+C*H-B*I-(C*E) +B*F) 
+AL2*(F*G-D*I-C*G+A*I+C*D-A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*B))/ 
(-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*D),3)) 


/((2*AL1*(POWER((-(F*H) + E*I),2)+ (C*H-B*I)*(F*G-D*I) + (-(C*E) +B*F)* 
(-(E*G) + D*H) + (-(F*H) + E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G)+A*I) + 

(-(C*E) + B*F)*(B*G-A*H) + (-(F*H) + E*I)*(-(C*E) +B*F) + (C*H-B*I)*(C*D-A*F) 
+(-(C*E)+B*F)*(-(B*D)+A*B)) 


+2*AL2*((F*G-D*I)*(-(F*H) +E*I) +(-(C*G) +A*I)*(F*G-D*I) + 
(-(E*G) +D*H)*(C*D-A*F) + (F*G-D*I)*(C*H-B*I) + 

POWER((-(C*G) +A*I),2)+(C*D-A*F)*(B*G-A*H) +(F*G-D*I)* 
(-(C*E) + B*F) + (-(C*G) + A*I)*(C*D-A*F) +(C*D-A*F)*(-(B*D)+A*B)) 
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+ 2*AL3*((-(E*G)+D*H)*(-(F*H) +E*D + (B*G-A*H)*(F*G-D*I) + 
(-(B*D) + A*E)*(-(E*G) +D*H) + (-(E*G) + D*H)*(C*H-B*I) + (B*G-A*H)* 
(-(C*G) + A*I) + (-(B*D) + A*E)*(B*G-A*H) + (-(E*G) + D*H)*(-(C*E) +B*F) 
+(B*G-A*H)*(C*D-A*F) +POWER((-(B*D)+A*E),2)))/ 


(POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*B*D),2)) 


-POWER(((AL1*(-(F*H) +E*I+C*H-B*I-(C*E) +B*F) 
+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 
+AL3*(-E*G +D*H+B*G-A*H-B*D+A*B))/ 
(-(C*E*G) + B*E*G + C*D*H - A*F*H - B*D*I + A*B*D),2))**1.5) - MM3),2); 
we 


* SUBJECT TO 
* 

EQ1.. 

ALI] =G= 0; 
EQ2?2.. 

AL2 =G= 0; 
EQ3.. 

AL3 =G= 0; 
EQ4.. 

ALI + AL2 + AL3 =L=1; 
EQS.. 

A =L=0>; 
EQ6.. 

E =L=0; 
EQ?7.. 

I =L=0; 
EQ8.. 

B =G=0; 
EQ9.. 

C =G=0; 
EQ10.. 

B+C =L=-A; 

EQ11 

D =G=0; 
EQ12 

F =G=0; 
EQ13.. 

D+F =L= -E; 

EQ14.. 

G =G=0; 
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EQ17.. 
-(AL1 *(-(F*H) + E*I+ C*H-B*I-(C*E) + B*F) 
+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 
+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 

(-(C*E*G) + B*F*G+C*D*H-A*F*H-B*D*I+A*E*l) =E= EX ; 


peenrneern BOR --------------------- aoe teed en 


* 


MODEL MATCHM /ALL/ ; 
* Parameters bounding 


A.LO = -7.0; B.LO = 3.0; C.LO = 0.0; D.LO = 0.0; E.LO =-7.0; 
A.UP = -5.0; B.UP = 5.0; C.UP = 2.0; D.UP = 1.0; E.UP =-5.0; 
F.LO = 4.0; G.LO = 0.0; H.LO = 0.0; I.LO = -8.0; 

F.UP = 5.0; G.UP = 1.0; H.UP = 2.0; I.UP = -7.0; 

AL1.LO = 1.0; AL2.LO = 0.0; AL3.LO = 0.0; 

AL1.UP = 1.0; AL2.UP = 0.0; AL3.UP = 0.0; 


* The initial solution 


A.L = AA; B.L = BB; C.L = CC; D.L = DD; EL = EE: 
F.L = FF;G.L=GG;HL=HH:IL=IQI; 
ALI.L = ALP1 ; AL2.L = ALP2 ; AL3.L = ALP3 ; 


SOLVE MATCHM USING NLP MINIMIZING Z ; 
DISPLAY Z.L,A.L , B.L, C.L, D.L, E.L, F.L, G.L, H.L, LL, 
ALI.L, AL2.L, AL3.L; 


* Write output to the CMS file 


FILE OUT! /FILE OBJ A/ ; 
PUT OUT] ; 

ror Z.L: 

FILE OUT2 /FILE ALPHA A/ ; 
PUT OUT?2 ; 

PUT A.L ; 

FILE OUT3 /FILE BRAVO A/; 
PUT OUTS3 ; 
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PUT B.L ; 

FILE OUT4 /FILE CHAREE A / ; 
PUT OUT4 ; 

PUT C.L ; 

FILE OUTS /FILE DELTA A/ ; 
PUT OUTS ; 

PUT D.L ; 

FILE OUT6 /FILE ECHO A/ ; 
PUT OUTS ; 

PUT E.L ; 


FILE OUT7 /FILE FOXTROT A / ; 


PUT OUT7 ; 

PUT EL ; 

FILE OUTS8 /FILE GOLF A/ ; 
PUT OUTS ; 

PUT G.L ; 

FILE OUT9Y /FILE HOTEL A / ; 
PUT OUT? ; 

PUR Hei 

FILE OUT10 /FILE INDIA A/ ; 
PUT OUTIO ; 

PUT L.L ; 

FILE OUT11 /FILE ALI A/; 
PUT OUTI1 ; 

PUT ALI1.L ; 

FILE OUT12 /FILE AL2 A/; 
PUT OUTI12 ; 

PUT AL2.L ; 

FILE OUT13 /FILE AL3 A/; 
PUT OUTI3 ; 

PUT AL3.L ; 
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THE INITIAL COMPLETE SYSTEM LIFETIME SET 1-5 


ow | sera | ser2 | ser3 | sers| ser s 
p2 | o.2455 | 0.2413 | 0.2257 | 0.1239 | 0.1999 

3 | 0.2662 
0.2925 
0.3205 
0.4214 
0.4394 
0.4685 0.3893 


0.5804 0.4761 0.4691 0.4686 0.4410 
0.5871 | 0.4909 | 0.5254 0.4766 


a0) 
ali 
a2 


0.6228 Oie5 2.910 0n55316 0.4885 0.4964 


ra 


0.6534 O58 8 05998 0.4970 0.5074 


a 
\9 


Oi. 7220 Oms9162 0.6575 0.5114 0.5240 
0.7590 8.6235 0.6623 06073 O02 5325 
0.7618 0.6419 OIN6855 0.6433 026033 


19 
16@ 
a8 
al 


0.8206 0.7156 | _0.7576 | 0.6255 
0.8245 0.7682 0.7690 0.8057 OF /39il 


~J 


0.8327 0.7723 
1.0056 0.9126 
1.0853 1.0832 
1.0853 1.1936 


nt 


N fe 
io |@ 


0 
1 
2 


NW Jd 


NO 
ee) 


NO 
He 


= 
a 


NO 
wn 


0.705 0.704 0.679 0.735 0.742 | 
0.70936 0.57558 0.80566 | 0.77290 
imiseenie) oPpaigaic i} 200991985) 14 40972°% ' 
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0.47917 
O'sGeai910 





THE INITIAL COMPLETE SYSTEM LIFETIME SET 6-10 


Tx [sere | ser? | sere] sere] oerao 
. ) 









0.6030 0.4523 
0.6281 0.5584 
il : 
ak 


12 


a 


Oi Sa6u 0.2775 0.3850 0.2977 0.4091 


0.6764 0.5800 
0.6851 0.5805 
0.7738 0.5993 
0.8541 0.6354 


1 
1 
18 
1 
1 
1 
18 
0 


WW 


0.8654 0.6261 0.3156 0.4055 0.6572 
0.4589 Orn6is 316 0.7584 0.4076 0.7529 
1098 0.6776 0.7677 0.4333 0.7874 


ce 


im0s 52 0.8705 0.8754 0.6178 0. 7991 
1.0679 1.1124 0.6836 0.6245 0.8572 


“J 


© 


1.1196 0.8920 
1.1995 | 1.1742 0.9024 
1.2303 0.9699 
1.3043 0.9734 
2.0273 1.0452 
3.2993 1.1463 
3.5360 1.1594 
1.000 0. 648 
0.82457 0.46844 | 


bt [dD 
ra 


2 
2 
20 
Z 


W |N 


U1 


2u™ | 


1 seorma4 1225687 | oOT79r02a | “reser7:t OL yaro ee 


| 


3° MM 


THE INITIAL COMPLETE SYSTEM LIFETIME SET 11-15 


Ver in | seraz | seras | seria | seri 
0.0268 
0.0745 | 0.1285 
0.1235 0.2052 0.1709 
0.1800 0.2435 0.1998 
0.2526 | 0.2924 0.2075 
0.2704 0.3335 0.2588 






es 


0.2972 | 0.3645 0.2669 

0.3697 | 0.4936 0.2737 

0.3811 0.4984 0.2508 

0.3990 0.5064 0.3506 

0.4070 | 0.5778 
0.4855 | 0.6699 
0.6402 | 0.7442 
0.7087 | 0.7502 
0.7408 | 0.8656 
0.8496 | 0.8831 
0.9290 | 1.0607 
0.9873 | 1.1003 
0.9983 | 1.1834 
1.0318 | 1.2235 0.7936 
1.0937 | 1.3370 1.1532 
Lgl 7 1.5754 1.2444 
1.2706 | 1.8085 1.2781 


1 
1 
2 
3 


he 1O 


3 


i 
1 
UL 
Aa 


@ 


2 


N te [eR 
Oo jw 


ND H 
=| 


NO 


al 
2 


bw Th 
Ww 


1.4875 | 1.9023 1.3428 
1.9894 1.5937 
0.665 0.831 | 9.560 | 0.872 | 0 46a 
0.64986 | 0.66778 0.68553 


1 0.32858 0.60712 | 0.69085} 1.51250} 0.77568 


NO Fd 


1.4958 


d MM 
3" 


THE INITIAL COMPLETE SYSTEM LIFETIME SET 16-20 


a 
1 | 0.1677 | 0.0565] 0.0398 | 0.1798 | 0.0202 | 
2 | 0.1819 | 0.1812] 0.1402 | 0.2037 | 0.1920 | 
he | 9.4393 | 14393 | 9.4147 ~4147 9.3787 oS. 8W7 | 9.3629 | S629 | 0.3645 | .3645 
ee 
=a 










N |e 


| 17 | 0.0390 | 0.7970] 0.8065 | 0.8492 | 0.7092 _| 
0.7447 


0 d ; : , 
ama | 0.ss- | ovcse | oom: | peel a 
2° 
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APPENDIX B. 
PROGRAM COST 


THIS PROGRAM PROVIDES THE OPTIMAL AGE REPLACEMENT OF THE 
PHASE TYPE DISTRIBUTION. IT WORKS ONLY SPECIFIC CASE OF 3 
PHASES OF TRANSIENT STATES. THE PARAMETERS ARE AS FOLLOWS: 
Td,J) = ENTRIES OF THE TRANSITION MATRIX OF PH DISTN. 
AL) = INITIAL PROBABILITY OF BEING IN STATE I 


Cl = PLANNED MAINTENANCE COST 
C2 = UNPLANNED MAINTENANCE COST 
16 = MAXIMUM LIFETIME THAT WANT TO CALCULATE 


DET = DETERMINANT OF TRANSITION MATRIX 
SMALL = THE OPTIMAL AGE REPLACEMENT 
CS = AGE REPLACEMENT COSTS 


INITIAL INPUT 


INTEGER I,J,K,N 
REAL C1, C2, TI(3,3), T(3,3), DET, AL(3), L, OBJ 
REAL TS(3,3), SMALL, PI, A, IDEN(3,3), S, NUMER 
REAL DENOM, P(3,3), Q(3,3), CS(400), EXPTS(3,3) 
DATA C1,C2,L,TS/100.0,500.0,4.0,9*0.0 / 

N =0 

S = 0.01 

OPEN(UNIT =11,FILE = ’ALPHA’) 

OPEN(UNIT =12,FILE = BRAVO’) 

OPEN(UNIT =13,FILE = ’CHAREE’) 

OPEN(UNIT =14,FILE = ’DELTA’) 

OPEN(UNIT =15,FILE = ’ECHO’) 

OPEN(UNIT =16,FILE = ’FOXTROT’) 
OPEN(UNIT =17,FILE = ’GOLF’) 

OPEN(UNIT =18,FILE = ’HOTEL’) 

OPEN(UNIT =19,FILE = ’INDIA’) 

OPEN(UNIT =20,FILE = ’AL1’) 

OPEN(UNIT =21,FILE = ’AL2’) 

OPEN(UNIT =22,FILE = ’AL3’) 

OPEN(UNIT =23,FILE = ’OBJ’) 

READ(11,*)T(1,1) 

READ(12,*)T(1,2) 

READ(13,*)T(1,3) 
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READ(14,*)T(2,1) 
READ(15,*)T(2,2) 
READ(16,*)T(2,3) 
READ (17,*)T(3,1) 
READ (18,*)T(3,2) 
READ(19,*)T(3,3) 
READ(20,*)AL(1) 
READ (21,*)AL(2) 
READ(22,*)AL(3) 
READ(23,*)OBJ 
IF(OBJ .GT. 0.05) GOTO 555 


* 


* CALCULATION COST FOR EACH REPLACEMENT TIME 
DET = -(T(1,3)*T(2,2)*T(3,1)) 

+ T(1,2)*T(2,3)*T(3,1) 

+ T(1,3)*T(2,1)*T(3,2) 

- T(1,1)*T(2,3)*T(3,2) 

- T(1,2)*T(2,1)*T(3,3) 

+ T(1,1)*T(2,2)*T(3,3) 
PRINT*,’DET = ’,DET 
TI(1,1) = (-(T(2,3)*T(3,2)) + T(2,2)*T(3,3))/DET 
T1(1,2) = (T(1,3)*T(3,2) - T(1,2)*T(3,3))/DET 
TI(1,3) = (-(T(1,3)*T(2,2)) + T1,2)*T(2,3))/DET 
T1(2,1) = (T(2,3)*T(3, 1) - T(2,1)*T(3,3))/DET 
T1(2,2) = (-(T(1,3)*T(3,1)) + T(1,1)*T(3,3))/DET 
T1(2,3) = (T(1,3)*T(2,1) - T(1,1)*T(2,3))/DET 
T1(3,1) = (-(T(2,2)*T(3,1)) + T(2,1)*T(3,2))/DET 
T1(3,2) = (T(1,2)*T(3,1) - T(1,1)*T(3,2))/DET 
T1(3,3) = (-(T(1,2)*T(2,1)) + T(1,1)*T(2,2))/DET 
DO 101 = 1,3 


— = “sa ™., 
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N = N+1 
DO 201 = 1,3 
DO 15J = 1,3 
TSC,J) = S*TC,J) 
15 CONTINUE 
20 CONTINUE 


%* 


* FIND THE EXPONENTIAL OF MATRIX TS 


*« 


CALL EXPON(TS,EXPTS, PI) 


A =0.0 
DO 211 = 1,3 
DO 22 J = 1,3 


A = A + AL()*EXPTS(,J) 
22 CONTINUE 
21 CONTINUE 
NUMER = C2 - A*(C2-C1) 
DO 301 = 1,3 
DO 25J = 1,3 
P(,J) = EXPTS(,J) - IDEN,J) 
25 CONTINUE 
30 CONTINUE 
DO 451 = 1,3 
DO 40J = 1,3 
Qd,J) = 0.0 
DO 35 K = 1,3 
Qd,J) = Qd,J) +TId,K)*P(K,J) 
35 CONTINUE 
40 CONTINUE 
45 CONTINUE 
DENOM = 0.0 
DO 551 = 1,3 
DO 50J = 1,3 
DENOM = DENOM + AL()*Q(,J) 
50 CONTINUE 
55 CONTINUE 
CS(N) = NUMER/DENOM 
S = $+0.01 
IF(S.LE.L)GOTO 100 


*« 


* FIND THE TIME THAT HAS MINIMUM COST 


*« 
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K =] 
200 IF(K.EQ.N)THEN 

SMALL = K*0.01 

ELSEIF(CS(K + 1).GT.CS(K)) THEN 
SMALL = K*0.01 

ELSEIF(CS(K + 1).LE.CS(K)) THEN 
SMALL = (K+1)*0.01 
K=K+1 
GOTO 200 

ENDIF 


*K 


* OUTPUT OPTIMAL REPLACEMENT AGE 
*« 
OPEN(UNIT = 3,FILE = ’OPTAGP’) 
WRITE(3,333)SMALL , PI 
333 FORMAT(1X,F6.3,4X,E11.4) 
OPEN(UNIT = 4,FILE = ’OPTCOST’) 
WRITE(4,*)CS(K) 
555 STOP 
END 


DK IK 2K 24 25 ac 4c 2 2c 2K 2 ic 2k aie 2k ic 2c 24 ik 24 ic 2K 24¢ 2k 28k 2c ic 28 ic 2k 2c ic akc 2c 2k 28 ic aK aC ac 2K 2iC 2K 28 2c aK 24C aK 2K 246 2c 2 ik ic 2k 24¢ 2K 24C 2K 2K 24¢ 2k 2 25k 2K 2c 2K 2k 2c 2K 2K 2k 


SUBROUTINE EXPON(A, EXPA, PI) 


THIS SUBROUTINE USES TO FIND THE EXPONENTIAL OF A 3X3 MATRIX. 
IT WORK WITH COUPLE OF IMSL SUBROUTINES 
CEVCRG’,’LINCG’,’EPIRG’). 


O2010°70' © 


INTEGER LDA,LDEVEC,N,NOUT,1LJ,K 
PARAMETER(N =3,LDA =N,LDEVEC=N,LDUINV =N,LD=N,LDU=N) 
REAL A(LDA,N),PI,EXPA(3,3) 
COMPLEX EVAL(N), EVEC(LDEVEC,N), ELD(3,3), UINV(3,3), 
COMPLEX C(3,3), D(3,3) 
OPEN(UNIT =1,FILE=’EXPONENT’) 
*K 
* CALCULATE THE EIGENVALUES AND EIGENVECTERS 
* 
CALL EVCRG(N,A,LDA,EVAL,EVEC,LDEVEC) 
DO 15I = 1,3 
DO 10J = 1,3 
ELDG,J = 0.0 
10 CONTINUE 
15 CONTINUE 
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ELD(1,1) = CEXP(EVAL(1)) 
ELD(2,2) = CEXP(EVAL(2)) 
ELD(3,3) = CEXP(EVAL(3)) 


ae 


* CALCULATE THE INVERSE EIGENVECTERS’ MATRIX 
* 
CALL LINCG(N,EVEC,LDU,UINV,LDUINV) 
DO 301 = 1,3 
DO 253 = 1,3 
Cd,) = 0.0 
DO 20 K = 1,3 
CdJ) = Cd,J) + ELD,K)*UINV(K,J) 
20 CONTINUE 
25 CONTINUE 
30 CONTINUE 
DO 451 = 1,3 
DO 40] = 1,3 
Dd,J) = 0.0 
DO 35 K = 1,3 
DJ) = Dd,J) + EVEC,K)*C(K,J) 
35 CONTINUE 
40 CONTINUE 
45 CONTINUE 
DO 551 = 1,3 
DO 50J = 1,3 
EXPA(I,J) = REAL(D(,J)) 
50 CONTINUE 
55 CONTINUE 
e 
C CALCULATE THE PERFORMANCE INDEX OF EIGENVALUES AND 
C EIGENVECTERS IF LESS THAN 1 ’EXCELLENT’, IF IT IS BETWEEN 1 
C AND 100 ’GOOD’, IF IT IS GREATER THAN 100 ’THE SOLUTION IS 
C SUSPECTED’ 
ec 
PI = EPIRG(N,N,A,LDA,EVAL,EVEC,LDEVEC) 
RETURN 
END 
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FQAAAANANINAINANAINIANM 


PROGRAM SIMLIF 


THIS PROGRAM PROVIDES A RANDOM NUMBER OF A PHASE TYPE 
DISTRIBUTION THAT HAVE SPECIFIC REPRESENTATION WITH AN 
INITIAL PROBABILITIES VECTOR(ALPHA) AND A MATRIX OF 
TRANSITION RATES BETWEEN TRANSIENT STATES. THIS PROGRAMS 
WORKS WITH SUBROUTINE ’RANNUM”’ AND MORE SPECIFIC CASE 
ONLY 4 STATES DISTRIBUTION. THE VARIABLES ARE AS FOLLOWS: 


LD(,J)= TRANSITION RATE FROM STATE I TO STATE J 

Pd,J) = TRANSITION PROBABILITY FROM STATE I TO STATE J 
AL() = INITIAL PROBABILITY OF BEING IN STATE I 

ELIF = THE SUPPOSED UPPER BOUND EXPECTED LIFETIME 


* INITIALIZATION 
xk 


*K 


INTEGER S,SEED,N,N1 
REAL U,EXP,LIF,ELIF 

REAL LD12,LD13,LD14,LD21,LD23,LD24,LD31,LD32,LD34 
REAL P12,P13,P21,P23,P31,P32,AL1,AL2,AL3,D 

DATA LD12,LD13,LD14,LD21 ,LD23,LD24,LD31,LD32,LD34 
+  /5.0,5.0,5.0,6.0,6.0,6.0,7.0,7.0,7.0/ 

DATA P12,P13,P21,P23,P31,P32,AL1,AL2,AL3,ELIF 

+  /0.8,0.2,0.167,0.833,0.143,0.286,1.0,0.0,0.0,10.0/ 

CALL EXCMS(’FILEDEF 10 DISK FILE LFTEST (DISP MOD’) 
OPEN(UNIT = 1,FILE =’LIFETIMEP’) 

OPEN(UNIT = 2,FILE =’SEED’) 

READ(2,*)SEED 

CLOSE(2) 

LIF = 0.0 


* GENERATE UNIFORM 0 AND 1 FOR INITIAL PROBABILITY 


* 


*K 


CALL RANNUM(1,SEED,0.0,1.0,0,U) 
IF (U.LE.AL1) THEN 
S=] 
ELSEIF(U.LE. AL] + AL2)THEN 
S=2 
ELSE 
S =3 
ENDIF 


* SIMULATE LIFETIME WHEN IS IN STATE 1 


* 
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5 CONTINUE 
IF ((S.EQ.1) .AND. (LIF.LT.ELIF)) THEN 
CALL RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P12+P13)THEN 
S=4 
CALL RANNUM(3,SEED,LD14,0,0,EXP) 
LIF = LIF + EXP 
ELSEIF(U.GT.P12)THEN 
S =3 
CALL RANNUM(3,SEED,LD13,0,0,EXP) 


Ss = 
CALL RANNUM(3,SEED,LD12,0,0,EXP) 
LIF = LIF + EXP 
ENDIF 
ENDIF 


* 


* SIMULATE LIFETIME WHEN IS IN STATE 2 
* 
IF((S.EQ.2) .AND. (LIF.LT.ELIF))THEN 
CALL RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P21+P23)THEN 
S=4 
CALL RANNUM(3,SEED,LD24,0,0,EXP) 
LIF = LIF + EXP 

ELSEIF(U.GT.P21)THEN 
S =3 
CALL RANNUM(3,SEED,LD23,0,0,EXP) 
LIF = LIF + EXP 

ELSE 
S=1 
CALL RANNUM(3,SEED,LD21,0,0,EXP) 
LIF = LIF + EXP 

ENDIF 

ENDIF 


* 


* SIMULATION LIFETIME WHEN IS IN STATE 3 
* 

IF((S.EQ.3) .AND. (LIF.LT.ELIF)) THEN 
CALL RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P31 +P32)THEN 

S=4 
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* 


CALL RANNUM(3,SEED,LD34,0,0,EXP) 
LIF = LIF + EXP 
ELSEIF(U.GT.P31)THEN 
S=2 
CALL RANNUM(3,SEED,LD32,0,0,EXP) 
LIF = LIF + EXP 
ELSE 
S=1 
CALL RANNUM(3,SEED,LD31,0,0,EXP) 
LIF = LIF + EXP 


IF((S.LT.4) .AND. (LIF.LT.ELIF))GOTO 5 


* OUTPUT TO ’FILE LIFETIME A’ 


WRITE(10,*)LIF 
WRITE(1,*)LIF 
OPEN(UNIT = 2,FILE =’SEED’) 
WRITE(2,*)SEED 

STOP 

END 


BK HE DHE DK DC DAC Die DIK DIC aie ae DIK DIC 24 aie DIC DIC dhe ie De DHE die BH DiC DIC de DC DiC DIC BIC DiC DIC DIK BHC DIC KC DC DHE DHE IHC dC DC BHC IHC DC DE DC IHC DIC DiC DC BK I4C D4 DiC IK IK IC IHC iC IK Dk D4 iC DiC IC dhe Dic dik IK IC 3kc 


QEQAAQQANNANNAaAIAANANM 


SUBROUTINE RANNUM(DISTN,SEED,RPARM1,RPARM2,IPARM,X) 


THIS SUBROUTINE PROVIDES A RANDOM NUMBER OF 7 DISTRIBUTION. 

IT INTERFACES WITH THE LLRANDOMII ROUTINES PROVIDED IN THE 

NONIMSL LIBRARY. THE PARAMETER AND CALLING PROCEDURE ARE 

AS FOLLOWS: 

DIST = DISTRIBUTION TYPE YOU WANT TO SELECT AN INTEGER 
BETWEEN 1 AND 7. 

SEED = THE RANDOM NUMBER SEED YOU WISH TO USE. 

RPARM1, RPARM2, AND IPARM ARE REAL AND INTEGER PARAMETERS. 

PASSED TO THE ROUTINE WITH MEANINGS WHICH VARY WITH THE 

TYPE OF DISTRIBUTION YOU DESIRE. 

X = THE RETURNED RANDOM NUMBER, IT IS ALWAYS REAL. 

DISTRIBUTION NUMBERS AND THE ASSOCIATED PARM DEFINITIONS 

1 = UNIFORM ON THE INTERVAL RPARM1 TO RPARM2 

2 = NORMAL WITH MEAN RPARMI AND VARIANCE RPARM2 

3 = EXPONENTIAL WITH RATE RPARM1 

4 = COUCHY WITH A = RPARMI AND B = RPARM2 

5 = GAMMA WITH SHAPE RPARM2 AND RATE RPARM1 
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C 6 = POISSON WITH RATE RPARM1 
C 7 = GEOMETRIC WITH P = RPARMI 
c 


REAL RPARMI,RPARM2,X 

INTEGER DISTN,SEED,IPARM,N 

REAL TEMP,VARIAT(1) 

IF(DISTN.LE.0 .OR. DISTN.GE.8)THEN 
WRITE(10,*)’DISTN DOES NOT EXIT” 
STOP 

ENDIF 


GOTO (10,20,30,40,50,60,70),DISTN 
GENERATE A UNIFORM BETWEEN RPARM1 AND RPARM2 


10 CONTINUE 
IF(RPARM1-RPARM2.EQ.0)THEN 
WRITE(10,*)’ ILLEGAL EQUAL RPARMS IN RANNUM’ 
STOP 
ENDIF 
IF(RPARM1.GT.RPARM2)THEN 
TEMP = RPARMI 


CALL LRND(SEED,VARIAT, 1,1,0) 
VARIAT(1) = RPARMI1 + (RPARM2-RPARM1)*VARIAT(1) 
GOTO 99 


GENERATE A NORMAL WITH MEAN RPARMI AND STD. DEV. RPARM2 


20 CALL LNORM(SEED, VARIAT, 1, 1,0) 
WRITE(10,*)NORMAL(0,1) ’, VARIAT(1) 
VARIAT(1) = (VARIAT(1)*RPARM2) + RPARMI1 
GOTO 99 


% 


GENERATE AN EXPONENTIAL WITH RATE RPARMI1 (1/MEAN) 


30 CONTINUE 
IF(RPARM1.EQ.0.0)THEN 
WRITE(10,*)" ILLEGAL ZERO RATE IN RANNUM? 
STOP 
ENDIF 
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% 


% 


CALL LEXPN(SEED,VARIAT, 1,1,0) 
VARIAT(1) = VARIAT(1)/RPARM1 
GOTO 99 


GENERATE A COUCHY WITH A = RPARM1 AND B = RPARM2 


40 CONTINUE 

IF(RPARM2.LE.0.0)THEN 
WRITE(10,*) ILLEGAL COUCHY SPREAD IN RANNUM, B= ’,RPARM2 
STOP 

ENDIF 

CALL LCCHY(SEED,VARIAT, 1,1,0) 

VARIAT(1) = (VARIAT(1)*RPARM2) + RPARMI 

GOTO 99 


GENERATE A GAMMA WITH SHAPE RPARM2 AND RATE RPARM1 


50 CONTINUE 

IF(RPARM1.LE.0.0) THEN 
WRITE(10,*)’ ILLEGAL NONPOSITIVE GAMMA RATE IN RANNUM’ 
STOP 

ENDIF 

IF(RPARM2.LE.0.0)THEN 
WRITE(10,*) ILLEGAL SHAPE PARAMETER IN RANNUM’ 
STOP 

ENDIF 

CALL LGAMA(SEED, VARIAT, 1,1,0,RPARMZ2) 

VARIAT(1) = VARIAT(1)*(1.0/RPARM1) 

GOTO 99 


GENERATE A POISSON WITH RATE RPARM1 


60 CONTINUE 
IF(RPARM1.LE.0.0)THEN 
WRITE(10,*) ILLEGAL POISSON RATE IN RANNUM’ 
STOP 
ENDIF 
CALL LPOIS(SEED,VARIAT,1,1,0,RPARM1) 
GOTO 99 


GENERATE A GEOMETRIC WITH P = RPARM1 


70 CONTINUE 
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IF(RPARM1.LE.0.0) THEN 
WRITE(10, *) ILLEGAL GEOM. PROB. IN RANNUM’ 
STOP 

ENDIF 

CALL LGEOM(SEED,VARIAT, 1,1,0,RPARM1) 

GOTO 99 


99 CONTINUE 
X = VARIAT(1) 
RETURN 
END 
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PROGRAM MOMENT 


THIS PROGRAM USES TO COMPUTE THE SECOND AND THIRD 
STANDARDIZED MOMENTS THAT SUPPORTS MOMENT MATCHING IN 
GAMS PROGRAM. 


e Ae) ©) C2 


ee DATA 


INTEGER LIMIT,N,N1,I,J,COUNT 
REAL LIF, NEW(2), REPAGE, A, B, C, EX, EXSQ, EXC, EX3 
REAL STD, SECOND, THIRD, C1, C2, AAC, OBJ 
REAL RAWDAT(1000,2), FB(1000), P(1000), TOT, TCOST 
DATA C1,C2,TOT, TCOST,AAC/100.0,500.0,0.0,0.0,0.0/ 
OPEN(UNIT = 1,FILE =’LIFETIME’ STATUS =’OLD’) 
OPEN(UNIT = 2,FILE =’OPTAGE’ STATUS =’OLD’) 
OPEN(UNIT = 3,FILE =’RAWDATA’ STATUS =’OLD’) 
OPEN(UNIT = 4,FILE =’MOMENT’) 
OPEN(UNIT = 5,FILE =’COUNT’,STATUS =’OLD’) 
READ(1,*)LIF 
READ(2,*)REPAGE 
READ(5,*)COUNT 
PRINT*,’COUNT=’ ,COUNT 
= 25 + COUNT 
DO 51 = 1,1000 
Pd) = 0.0 
FB() = 0.0 
5 CONTINUE 
DO 441 = 1,N1-1 
READ(3,11)RAWDAT(,1),RAWDAT(,2) 
11 FORMAT(1X,F7.4,4X,F3.1) 
44 CONTINUE 
CLOSE(3) 
OPEN(UNIT =3,FILE=’RAWDATA’) 
we 


* CALCULATION 


* 


IF(LIF.LE.REPAGE)THEN 


NEW(1) = 

NEW(2) = 1.0 
ELSE 

NEW(1) = REPAGE 

NEW(2) = 0.0 
ENDIF 
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CALL INSERT(N1 ,NEW,RAWDAT) 
DO 55I = 1,NI 
WRITE(3,33) RAWDAT(, 1), RAWDAT(,2) 
33 FORMAT(1X,F7.4,4X,F3.1) 
55 CONTINUE 
IF(RAWDAT(1,2) .NE. 0.0)THEN 
FB(1) = (N1-1.0)/N1 


DO 10I = 2,NI 
IF(RAWDAT(L,2).EQ.0.0) THEN 
FB() = FB(I-1) 
ELSEIF(I.LT.N1)THEN 
FB(I) = ((N1-D/(N1 +1.0-1))*FB(I-1) 
ELSE 


P(1) = 1.0-FB(1) 
DO 22 I = 2,NI1 
P() = FB(-1)-FB() 
22 CONTINUE 
EX = 0.0 
EXSQ = 0.0 
EXC = 0.0 
EX3 = 0.0 
DO 15J = 1,N1 
TOT = TOT + RAWDAT(J,1) 
TCOST = TCOST + C2*RAWDAT(J,2) + C1*(1-RAWDAT(J,2)) 
A = RAWDAT(,1)*P(J) 
B = (RAWDAT(,1)**2)*P(J) 
C = (RAWDAT(J,1)**3)*P(J) 
EX = EX +A 
EXSQ = EXSQ + B 
EXC = EXC +C 
15 CONTINUE 
AAC = TCOST/TOT 
STD = (EXSQ - EX**2)**0.5 
EX3 = EXC + (2*EX**3) - (3*EX*EXSQ) 
SECOND = STD/EX 
THIRD = EX3/(STD**3) 
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*K 


* OUTPUT 
* 
WRITE(4, 100)SECOND, THIRD, AAC, EX 
100 FORMAT(2X,F11.8,4X,F11.8,4X,F8.3,4X,F6.3) 
*  CLOSE(5) 
*  QPEN(UNIT = 5,FILE =’COUNT’,STATUS =’OLD’) 
*  WRITE(5,*)COUNT+1 
STOP 
END 


DK 3K DIK IK SC IC SIC IC DIK 3K IC SC IIE DIC 3K DIC DAC IK 24C DIC IC DK BHC IC Dk DiC DIK DAC iC DAC DIK IC DE DE DIC 4c ic 2c DHE DHE afc DIC ic DAC DIC IC IC 3K IC 24K DIC iC DHE 2kC IC DIC iC DIK DIC kc IC 24 2c 2K 24 2c 2c 2k 


SUBROUTINE INSERT(LIMIT,NEW,X) 
C THIS SUBROUTINE USES TO INSERT A PAIR ’NEW’ DATA INTO THE 
C DATA FILE ’X’ IN ASCENDING ORDER. 

INTEGER LIMIT,J,K 

REAL NEW(2),X(1000,2) 

LOGICAL DONE 

DONE = .FALSE. 

1 — 

5 IF(J.LT.LIMIT).AND.(.NOT.DONE))THEN 
IF(X(J, 1). LT. NEW(1)) THEN 
J =J+1 


IF(J.EQ.LIMIT)THEN 
X(J,1) = NEW(1) 
X(J,2) = NEW(2) 
ELSE 
IF(J. LT. LIMIT)THEN 
DO 10 K = LIMIT,J+1,-1 
X(K, 1) = X(K-1,1) 
X(K,2) = X(K-1,2) 
10 CONTINUE 
X(J,1) = NEW(1) 
X(J,2) = NEW(2) 
ENDIF 
ENDIF 
RETURN 
END 
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<< *< Oee@remne Gl cme CY (2 2 12 0) @ 


PROGRAM DATA 


THIS PROGRAM PROVIDES THE RESULT SEQUENTIAL ESTIMATED 
ENTRIES TRANSITION RATES T(,J) AND INITIAL STATE PROB. 
ALPAH(),AND ALSO ADJUST, PREPARE ALL THE INITIAL DATA THAT 
WILL BE USED FOR THE NEXT SEQUENTIAL ESTIMATION. 

THE VARIABLE ARE AS FOLLOWS: 

A,B,C 

D,E,F = THE ESTIMATED TRANSITION MATRIX 

G,H,I 

AL1,AL2,AL3 = THE ESTIMATED INITIAL STATES PROBABILITIES 
X = RANDOM VARIABLE OF LIFETIME 

OPTAGE = OPTIMAL AGE REPLACEMENT 

SECOND = THE SECOND STANDARDIZED MOMENT 

THIRD = THE THIRD STANDARDIZED MOMENT 

OBJ = OBJECTIVE FUNCTION VALUE OF GAMS PROGRAM 


INPUT 


INTEGER COUNT 

REAL A,B,C,D,E,F,G,H,1,AL1,AL2, AL3,X 

REAL EX,OPTCS,OPTAGE,Z,SECOND, THIRD, DI,OBJ,AAC,PI 
OPEN(UNIT = 10,FILE =’OBJ’,STATUS =’OLD’) 

OPEN(UNIT = 11,FILE =’ALPHA’,STATUS =’OLD’) 
OPEN(UNIT = 12,FILE =’BRAVO’, STATUS =’OLD’) 
OPEN(UNIT = 13,FILE =’CHAREE’ STATUS =’OLD’) 
OPEN(UNIT = 14,FILE =’DELTA’ STATUS =’OLD’) 
OPEN(UNIT = 15,FILE =’ECHO’ STATUS =’OLD’) 
OPEN(UNIT = 16,FILE =’FOXTROT’,STATUS =’OLD’) 
OPEN(UNIT = 17,FILE =’GOLF’ STATUS =’OLD’) 
OPEN(UNIT = 18,FILE =’HOTEL’,STATUS =’OLD’) 
OPEN(UNIT = 19,FILE =’INDIA’,STATUS =’OLD’) 
OPEN(UNIT = 20,FILE =’AL1’,STATUS =’OLD’) 

OPEN(UNIT = 21,FILE =’AL2’,STATUS =’OLD’) 

OPEN(UNIT = 22,FILE =’AL3’,STATUS =’OLD’) 

OPEN(UNIT = 23,FILE =’OPTAGE’,STATUS =’OLD’) 
OPEN(UNIT = 24,FILE =’LIFETIME’ STATUS =’OLD’) 
OPEN(UNIT = 25,FILE =’MOMENT’,STATUS =’OLD’) 
OPEN(UNIT = 26,FILE =’COUNT’,STATUS =’OLD’) 
OPEN(UNIT = 27,FILE =’OPTCOST’ STATUS =’OLD’) 

CALL EXCMS(’FILEDEF 50 DISK FILE ESTDAT1 (DISP MOD’) 
CALL EXCMS(’FILEDEF 60 DISK FILE ESTDAT2 (DISP MOD’) 
CALL EXCMS(’FILEDEF 70 DISK FILE ESTDAT3 (DISP MOD’) 
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OPEN(UNIT = 30,FILE =’SCALAR2’) 
READ(10,*)OBJ 
IF(OBJ .GT. 0.05) THEN 

A = -1.0 


READ(11,*)A 
READ(12,*)B 
READ(13,*)C 
READ(14,*)D 


READ(20,*)ALI 

READ(21,*)AL2 

READ(22,*)AL3 
READ(23,*)OPTAGE,PI 

READ(24,*)X 
READ(25,*)SECOND, THIRD, AAC,EX 
READ(26,*)COUNT 
READ(27,*)OPTCS 

CLOSE(26) 

PRINT*,’OBJ=’,OBJ 


IF(OPTAGE .GT. X)THEN 
Z=xX 
DI = 1.0 
ELSE 
Z = OPTAGE 
DI = 0.0 
ENDIF 
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* 


* OUTPUT 
* 
OPEN(UNIT = 26,FILE =’COUNT’ STATUS =’OLD’) 
WRITE(26,*)COUNT+1 
WRITE(50,200)COUNT,A,B,C,D,E,F,G,H,I 
WRITE(60,300)COUNT,OPTCS, OPTAGE, X,Z,DI,AAC 
WRITE(70, 100)COUNT,AL1,AL2,AL3,SECOND, THIRD, OBJ, PI 


* 


* PREPARE THE DATA FOR "GAMS" PROGRAM 


* 


WRITE(30,*)’ SCALAR’ 


WRITE(30,400)’ AA INITIAL SOLUTION OF A /’,A,’ /” 
WRITE(30,400)’ BB INITIAL SOLUTION OF B /’,B,’ /” 
WRITE(30,400)’ CC INITIAL SOLUTION OF C /’,C,’ /’ 
WRITE(30,400)’ DD INITIAL SOLUTION OF D /’,D,’ /” 
WRITE(30,400)’ EE INITIAL SOLUTION OF A /’,E,’ /’ 
WRITE(30,400)’ FF INITIAL SOLUTION OF B/’,F,’ /” 
WRITE(30,400)’ GG INITIAL SOLUTION OF C /’,G,’ /’ 
WRITE(30,400)’ HH INITIAL SOLUTION OF D /’,H,’ /’ 
WRITE(30,400)’ II INITIAL SOLUTION OF D /’,I,’ /’ 
WRITE(30,500)’ ALP1 INITIAL SOLUTION OF ALI /’,AL1,’ /’ 
WRITE(30,500)’ ALP2 INITIAL SOLUTION OF AL2 /’,AL2,” /’ 
WRITE(30,500)’ ALP3 INITIAL SOLUTION OF AL3 /’,AL3,” /” 
WRITE(30,600)’ MM2 SECOND STANDARD MOMENT /’,SECOND,’/’ 
WRITE(30,600)’ MM3 THIRD STANDARD MOMENT /’ ,THIRD,’/’ 
WRITE(30,700)’ EX EXPECTED VALUE OF LIFETIME /’,EX,’ / ;’ 


100 FORMAT(1X,I3,3(2X,F5.2),2(2X,F8.5),2X,F4.2,2X,E11.4) 
200 FORMAT(1X,13,2X,9(2X,F6.2)) 
300 FORMAT(1X,13,2X,F8.3,2X,F5.2,2(4X,F7.4),4X,F3.1,4X,F7.3) 
400 FORMAT(A35,F6.2,A2) 
500 FORMAT(A39, F5.2,A2) 
600 FORMAT(A36,F11.8,A1) 
700 FORMAT(A40,F6.3,A4) 
STOP 
END 
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